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ABSTRACT 

We consider the problem of reconstructing time sequences of spa- 
tially sparse signals (with unknown and time-varying sparsity pat- 
terns) from a limited number of linear "incoherent" measurements, 
in real-time. The signals are sparse in some transform domain re- 
ferred to as the sparsity basis. For a single spatial signal, the solu- 
tion is provided by Compressed Sensing (CS). The question that we 
address is, for a sequence of sparse signals, can we do better than 
CS, if (a) the sparsity pattern of the signal's transform coefficients' 
vector changes slowly over time, and (b) a simple prior model on 
the temporal dynamics of its current non-zero elements is available. 
The overall idea of our solution is to use CS to estimate the sup- 
port set of the initial signal's transform vector. At future times, run 
a reduced order Kalman filter with the currently estimated support 
and estimate new additions to the support set by applying CS to the 
Kalman innovations or filtering error (whenever it is "large"). 

Index Terms/Keywords: compressed sensing, Kalman filter- 
ing, compressive sampling, sequential MMSE estimation 

1. INTRODUCTION 

We consider the problem of reconstructing time sequences of spa- 
tially sparse signals (with unknown and time-varying sparsity pat- 
terns) from a limited number of linear "incoherent" measurements, 
in real-time. The signals are sparse in some transform domain re- 
ferred to as the "sparsity basis" (I]. A common example of such a 
problem is dynamic MRI or CT to image deforming human organs 
such as the beating heart or to image brain neural activation patterns 
(in response to stimuli) using fMRI. The ability to perform real-time 
MRI capture and reconstruction can make interventional MR practi- 
cal ||2]. Human organ images are usually piecewise smooth and thus 
the wavelet transform is a valid sparsity basis |T][3l. Due to strong 
temporal dependencies, the sparsity pattern usually changes slowly 
over time. MRI captures a small (sub-Nyquist) number of Fourier 
transform coefficients of the image, which are known to be "inco- 
herent" with respect to the wavelet transform ||T][3]. Other examples 
include sequentially estimating optical flow of a single deforming 
object (sparse in Fourier domain) from a set of randomly spaced op- 
tical flow measurements (e.g. those at high intensity variation points 
l4l). or real-time video reconstruction using single-pixel camera l5l . 

The solution to the static version of the above problem is pro- 
vided by Compressed Sensing (CS) ||T||6l|2|. The noise-free ob- 
servations case IXI is exact with high probability (w.h.p.) while the 
noisy case Q has a small error w.h.p.. But existing solutions for the 
dynamic problem OH) treat the entire time sequence as a single spa- 
tiotemporal signal and perform CS to reconstruct it. This is a batch 
solution (need to wait to get the entire observation sequence) and has 
very high complexity. An alternative would be to apply CS at each 
time separately, which is online and low-complexity, but will require 
many more measurements to achieve low error (see Fig. [T). The 
question that we address is: can we do better than performing CS at 



each time separately, if (a) the sparsity pattern (support set) of the 
transform coefficients ' vector changes slowly, i. e. every time, none 
or only a few elements of the support change, and (b) a simple prior 
model on the temporal dynamics of its current non-zero elements is 
available (but do not know which coordinates are non-zero). 

Our solution is motivated by reformulating the above problem 
as causal minimum mean squared error (MMSE) estimation with a 
slow time-varying set of dominant basis directions (or equivalently 
the support of the transform vector). If the support is known, the 
MMSE solution is given by the Kalman filter (KF) |[9] for this sup- 
port. But what happens if the support is unknown and time- varying? 
The initial support can be estimated using CS (7). If at a given time, 
there is an addition to the support set, but we run the KF for the old 
model, there will be a model mismatch and the innovation (and fil- 
tering) error will increase. Whenever it does, the change in support 
can be estimated by running CS on the innovation or the filtering 
error, followed by thresholding. A Kalman update step is run using 
the new support set. If some coefficients become and remain nearly 
zero (or nearly constant), they can be removed from the support set. 

If, for a moment, we assume that CS Q gives the correct esti- 
mate of the support at all times, then the above approach will give 
the MMSE estimate of the signal at all times. The reason it is very 
likely that CS Q gives the correct estimate is because we use it to 
fit a very sparse "model change" signal to the filtering error. Also 
note that a full Kalman filter 15], that does not use the fact that the 
signal is sparse, is meaningless here, because the number of observa- 
tions available is smaller than the signal dimension, and thus many 
elements of the signal transform will be unobservable. Unless all 
unobservable modes are stable, the error will blow up. Other re- 
cent work that also attempts to use prior knowledge with CS, but to 
reconstruct only a single signal is II0IIIl|[T2l . 

2. THE MODEL AND PROBLEM FORMULATION 

Let {zt)mxi denote the spatial signal of interest at time t and 
{yt)nxi, with n < m, denote its observation vector at t. The 
signal, 2t, is sparse in a given sparsity basis (e.g. wavelet) with 
orthonormal basis matrix, "I>mxm, i.e. Xt — zt is a sparse vec- 
tor (only St « m elements of Xt are non-zero). Here ' denotes 
transpose. The observations are "incoherent" w.r.t. the sparsity basis 
of the signal, i.e. yt — Hzt + uit = H^xt + wt, where H„xm 
is such that the correlation between the columns of A = H$ is 
small enough to ensure that any St-column sub-matrix of A is "ap- 
proximately orthonormal" (its nonzero singular values are between 
to ^/TT5 for 5 < 1). Wt is the measurement noise. Thus 
the measurement model is: 

yt = Axt + wt, A^ H^, Wt ~ AA(0, a^^) (1) 

with Wt being temporally i.i.d.. We refer to xt as the state at t. Our 
goal is to get the "best" causal estimate of xt (or equivalently of 
zt = '^Xt) at each t. 



Let Tt denote the the support set of Xt, i.e. the set of its non- 
zero coordinates and let St = size{Tt). In other words, Tt = 
[ii, 12, . . . ist] where ik are the non-zero coordinates of Xt- For any 
set T, let {v)t denote the size{T) length sub-vector containing the 
elements of v corresponding to the indices in the set T. For another 
set, 7, we also use the notation T-, which treats T as a vector and 
selects the elements of T corresponding to the indices in the set 7. 
For a matrix A, At denotes the sub-matrix obtained by extracting 
the columns of A corresponding to the indices in T. We use the no- 
tation {Q)ti,t-2 to denote the sub-matrix of Q containing rows and 
columns corresponding to the entries in Ti and T2 respectively. The 
set operations U, H, and \ have the usual meanings (note Ti \ T2 
denotes elements of Ti not in T2). We use ' to denote transpose. T'^ 
denotes the complement of T w.r.t. [1 : m], i.e. T'^ = [1 : m] \ T. 
Also I |u| |p is the Ip norm of the vector v, i.e. \ \v\\p = \ vi\'')^^^' ■ 

Assumption 1. We assume slow changes in sparsity patterns, 
i.e. the maximum size of the change in the support set at any time 
is smaller (usually much smaller) than 5** at any t, i.e. Sm ff,max — 
iaa,xt[size{Tt \ Tt-i) + size{Tt-i \ Tt)] < mint St- This, as 
we shall see later, ensures that performing CS in Gaussian noise to 
only estimate the additions to the support set results in smaller error 
w.h.p. than doing CS to estimate the entire Tt again. 

Assumption 2. We also assume that A satisfies the Uniform 
Uncertainty Principle (UUP) (equation 1.6 of Q) at a sparsity level, 
Smax — maxt St- At first thought, it may appear that having A 
satisfy UUP at the smaller level Sdiff,max is sufficient, but, it is not. 
This is because if St > Sdiff.max for some t, even though Art is 
a tall matrix, there is no guarantee that at least St of its columns are 
linearly independent, i.e. it may happen that rank{ATt) < St- In 
this case, the system no longer remains observable for some elements 
of Xt- If any of these elements follow an unstable prior dynamic 
model, the KF error will increase unboundedly with t- 

System Model for xt- For the currently non-zero coefficients of 
Xt, we assume a spatially i.i.d. Gaussian random walk model, with 
noise variance aly^. At the first time instant at which {xt)i becomes 
non-zero, it is assumed to be generated from a zero mean Gaussian 
with variance (J^„it . Thus, we have the model: xo = 0, 

{xt), = {xt^i)i + {vt),, {vt)i ~ A/'(0,afy,), if i G Tf, i G T-i 
{xt), = (it-i), + {vt)., {vt)^'-^ N{Q,al,u) if i £Tt,ii T-i 
{Xt)r = {xt-i)^\fiiTt (2) 

The above model can be compactly written as: xq — 0, 

Xt = xt~i + i^t, i^t JV(0,Qt), 
(Qt)TtnTt_i,TtnTt_i = cr^ysl 

(Qt)Tt\Tt_i,Tt\Tt_i = <^'initl 

(Ot)Tf ,Tf = (3) 

where the set Tt is unknown Vi. If Tt were known at each t, i.e. 
the system model was completely defined, the MMSE estimate of 
Xt from yi,y2, . . .yt would be given by a reduced order KF defined 
for {xt)Tt- But, as explained in Sec. 1, in most practical problems, 
Tt is in fact unknown and time-varying. Often, it may be possible 
to get a rough prior estimate of Ti by thresholding the eigenvalues 
of the covariance of xi (possible to do if multiple realizations of 
xi are available to estimate its covariance). But without multiple 
i.i.d. realizations of the entire {xt}, which are impossible to obtain 
in most cases, it is not possible to get a-priori estimates of T for all 
t. But note that, it is possible to estimate <j1yg,a1„it for the model 
of l[3j using just one "training" realization of {xt} (which is usually 



easy to get) by setting the near-zero elements to zero in each xt and 
using the rest to obtain an ML estimate. 

Assuming known values of crly^, o-f^^f , our goal here is to get 
the best estimates of Tt and Xt at each t using y\, . . .yt- Specifically, 

1. At each time, t, get the best estimate of the support set, Tt, 
i.e. get an estimate Tt with smallest possible [size{Tt \ Ti) + 
size{ft \ Tt)] using yi,j/2 ---yt- 

2. Assuming the estimates of T\, .. .Tt are perfect (have zero 
error), get the MMSE estimate of Xt using y\,y2 . ■ .yt- 

3. KALMAN FILTERED COMPRESSED SENSING (KF-CS) 

We explain the proposed KF-CS algorithm below. 

Running the KF- Assume, for now, that the support set at the 
first time instant, Ti, is known. Consider the situation where the first 
change in the support occurs at a t = ta, i.e. for t < ta, Tt = 
Ti, and that the change is an addition to the support. This means 
that for t < ta, we need to just run a regular KF, which assumes 
the following reduced order measurement and system models: yt = 
AT{xt)T + wt, {xt)T = {xt-i)T + {vt)T, wlth T = Ti . The KF 
prediction step for this model is (9|: 

it\t-i = xt-i 

{Pt\t-l)T.T = {Pt-l)T,T + ^lysl (4) 

with initialization io = Q, Po = 0. The update step is (9): 

K ^ (Pt|t_i)T,T^TS-\ E.e ^ AT{Pt\t-l)T,TA'T+al,J 

(xt)T = (it|t.^i)T + K[yt - Axt\t-i] 

{xt)T'= = {Xt\t-l)T'= = {xt-ljT" 
{Pt)T,T =[I- KAT]{Pt\t-l)T,T (5) 

Note the update step for the elements of T'^ . 

Detecting If Addition to Support Set Occurred. The Kalman 
innovation error is yt — yt — Axt\t-i- For t < ta, yt = A{xt — 
a;i|t_i) + uit ~ A/'(0, Eie) (9). At t = ta, a new set. A, gets added 
to the support of a::t,i.e. yt = y4T(a;t)T + AA(2;t)A + wt, where the 
set A is unknown. Since the old model is used for the KF prediction, 
yt SLlt — ta will have non-zero mean, AA{xt)A, i.e. att = ta, 

yt = AA{xt)A + Wt, Wt = [Arixt - Xt\t-l)T + Wt] ~ A/'(0, E^e) (6) 

Thus, the problem of detecting if a new set has been added or not gets 
transformed into the problem of detecting if the Gaussian distributed 
yt has non-zero or zero mean. We use the well-known generalized 
Likelihood Ratio Test for this problem, which detects if the weighted 
innovation enw norm, I EN = y^Yl'^^^yt ^ threshold. One can 
replace lEN by the weighted norm of the filtering error, defined 
below, which will make the detection more sensitive. 

Estimating the Additions (using CS)- If the lEN (or FEN) is 
"high" (above the detection threshold), there is a need to estimate A. 
Att — ta, the innovation error, ijt, can also be written as: 

yt = AT'={xt)T'= +wt, Wt ~7V(0,E„) (7) 

while the filtering error, ijt.f — yt — Axt, can be written as: 

ytj = AT<^{xt)T'-- + Arixt — Xt)T + Wt 

— [I — ATK]AT<:{xt)T<= + Wt.f, WtJ =\I — ATK]wt 
WtJ ~ AA(0, Eye), S/e = [/ - ATA']E,e[/ - AtK]' (8) 



Algorithm 1 Kalman Filtered Compressive Sensing (KF-CS) 



Initialization: Set xo — 0, Po — 0, To= empty (if unknown) or equal to the known/partially known support. For t > 0, do, 

1. SetT ^ Tt-i. 

2. KF prediction. Run {4) using the current T. 

3. KF update. Run {Sjl using the cuiTent T. 

4. Addition (using CS). Compute lEN = y't'E~}yt where yt — yt — ATj|f_i (or compute FEN = y^j'T,J^yt.f where yt.f — 
yt — Axt), and check if it is greater than its threshold. If it is, 

(a) Run CS on the filtering error, yt.f — yt — Axt, i.e. compute the Dantzig selector using 

(b) Compute the support set of Xnew by thresholding, i.e. compute nz = {i : \{xnew)i\ > a}. Then the addition to the support set 
of is A = {T'^)nz- The new support set for xt is Tnew = T U A. 

(c) Set (Pt|t-i)A,A = crf„it/. Set T ^ T„ew 

(d) Run the KF update given in Q for the current T. 

Performance can be improved by iterating the above four steps until si2:e(A) = or FEN less than its threshold. 

5. Deletion. Compute the set A_d ~ {i &T : X]t=t-fc+i(^^)i ^ ka^}. The new support set is Tnew = T \ Ad- 

(a) Set (Pt|t_i)Aj3,[l:m] = 0, {Pt\t-l)[l:,n],Ao = . Ssl T ^ T„ew ■ 

6. KF update. Run {Sjl for the current T. 

7. Assign Tt ^ T. Output Tt, xt and the signal estimate, zt = $£4. Increment t and go to step[T] 



with (xt)T<: being a sparse vector with support, nz, s.t. A = 
(r'^)nz. Note that E/e < Eie. The problem of estimating {xt)T'^ 
from either yt or ytj is the problem of compressed sensing in 
Gaussian noise studied in (7J, except that the noise now is spatially 
colored. It is not immediately clear whether to use yt or ytj. In 
yt.f, the "noise", ilh.f, is smaller (the change {xt — xt-i)T has been 
estimated and subtracted out), but the new component, At" {xt)T'' , 
is also partially suppressed. But the suppression is small because 

ATKAT-{xt)T'' = AT{Pt]t_^alts + A'TAT)'^A'TAT-=ixt)T- 

(follows by rewriting K using matrix inversion lemma (9)) and 
A'j,AT<:{xt)T<: ~ A'rpAA{xt)A is small (since A satisfies UUP at 
level Smax and St^ = size{T U A) < Smax)- Thus, we use yt.f. 

The Dantzig selector |7] can be applied to yt.f to estimate 
(a;t)T<:, followed by using thresholding to compute its support, nz, 
and thus get A. Since wt.f is colored, one way to modify the 
Dantzig selector to apply it to l[8j is to compute: 

Xnew=SiTgmm\\f3\\i,s.t.\\A'TcA^^^''^U'[yt.f - At^P]\\oo < A™, (9) 

/3 

where E/e = UAU' is the eigenvalue decomposition of E/e and 
Am = -^2 log m trace (E/e)/?!. The set nz is estimated as nz — 
{i : \{xnew)i\ > ce} for some zeroing threshold a and the set A 
is A = {T'')„^. Thus the new support set is T„ew = T U A. We 
initialize the prediction covariance along A: (Pt|t_i)A,A = (j'i„itl. 

KF Update. We run the KF update given in (|5j with T = T„sm. 
This can be interpreted as a Bayesian version of Gauss-Dantzig (T). 

Iterating CS and KF-update. Often, it may happen that not all 
the elements of the true A get estimated in one run of the CS step, 
because of the error in CS. To address this, the above steps (CS and 
KF update) can be iterated until FEN goes below a threshold or un- 
til the estimated set A is empty or for a fixed number of iterations. A 
similar idea forms the basis of iterative CS reconstruction techniques 
such as stagewise Orthogonal Matching Pursuit 1 1 31 . 

Deleting Near-Zero Coefficients. Over time, some coefficients 
may become zero (or nearly-zero) and remain zero. Alternatively, 
some coefficients may wrongly get added, due to CS error. In both 



cases, the coefficients need to be removed from the support set Tt. 
A simple way to do this would be to check if \{xt)i\ < a (a is 
a zeroing threshold) for a few time instants. When a coefficient is 
removed, the corresponding row and column in Pt\t—i is set to zero. 

Deleting Constant Coefficients. If a coefficient, i, becomes 
nearly-constant (this may happen in certain applications), one can 
keep improving the estimate of its constant value by changing the 
prediction step for it to {Pt\t-i)i,i ~ {Pt-i)i,i. Either one can keep 
doing this forever (the error in its estimate will go to zero with t) 
or one can assume that the estimation error has become negligibly 
small after a finite time and then remove the coefficient index from 
Tt. The bias- variance tradeoff needs to be evaluated before select- 
ing either scheme. If the latter is done, then for future times, one 
needs to replace yt by yt - Ai{x)t and set (Pt|t-i)i,[i:m] = 0, 
{Pt\t-i)[i:m],i = 0. This will be part of future work. 

Initialization. Initially, the support set, Ti may be roughly 
known (estimated by thresholding the eigenvalues of the covariance 
of xi, which is computable if its multiple realizations are available) 
or unknown. We initialize KF-CS by setting xq = 0, Pq — and 
To = roughly known support or To = empty (if support is completely 
unknown). In the latter case, automatically alt — 1, the I EN (or 
FEN) will be large, and thus CS will run to estimate Ti . 

The entire KF-CS algorithm is summarized in Algorithm 1. 

3.1. Discussion 

The key difference between KF-CS and regular CS at each t is that 
KF-CS performs CS on the filtering error, yt.f, to only detect new 
additions while regular CS performs CS on the observation, yt, to 
detect the entire vector Xt (without using knowledge of the previ- 
ous support set). From (7] Theorem I.I] (the result should hold di- 
rectly or with slight modification for CS in colored noise, but we 
have not checked it yet), the CS error is directly proportional to the 
sparsity size of the vector being estimated and to the noise variance. 
But the dependence on sparsity size is much stronger (highly non- 
linear) while that on the noise variance is linear. Note that "CS on 
filtering error" needs to estimate at most Sdi f /,max coefficients and 
if f, max < Sti (Assumption 1). Thus, if the "noises" in both 
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Fig. 1. MSB plots of KF-CS (labeled CS-KF in the plots), Ti unknown and known cases, compared against regular CS in the first 3 figures 
and against the Full 256-dim KF in the last figure (its MSB is so large that we cannot plot it in the same scale as the others). The benchmark 
(MMSE estimate with known Ti, Ts) is the genie-aided KF. The simulated signal's energy at t is E[| |a;t| I2] ~ Siaf^^^ + (X]t=2 'S'r)o"sj,s). 



cases were the same, then the error in "CS on filtering error" would 
be much smaller than that in CS. But the "noise" in yt is only wt, 
while the "noise" in ijtj is wt plus Arixt — xt)T- If we can show 
stability of KF-CS error, and if the temporal prior is strong enough, 
the covariance of this extra "noise" will be small enough. Thus, it 
should be possible to show that the error in KF-CS is smaller than 
that in regular CS, if the sparsity pattern changes slowly enough and 
the temporal prior is strong enough. 

A rigorous comparison of KF-CS with CS and an analysis of 
KF-CS error stability will be part of future work. We provide here 
a qualitative discussion of the sources of error in KF-CS. The error 
can be due to two reasons - either the CS step misses some non- 
zero coefficients (or the deletion step wrongly deletes a non-zero 
coefficient) or the CS step estimates some extra coefficients (or the 
deletion step misses removing some coefficients). 

Occasionally, an element, i, of the true A may get missed, due 
to nonzero CS error (or the thresholding error). The extra error due 
to an {xt)i getting missed cannot be larger than the CS error at the 
current time (which itself is upper bounded by a small value w.h.p. 
171 ) plus a (due to thresholding). Also, eventually, when the mag- 
nitude of {xt)i increases at a future time, it will result in a "high" 
I EN and FEN, and the CS step at that time will estimate it w.h.p.. 

We can prevent too many extra coordinates from getting wrongly 
estimated by having a rough idea of the maximum sparsity of Xt and 
using thresholding to only select that many, or a few more, highest 
magnitude non-zero elements. The deletion threshold also needs to 
be selected appropriately. Also, if, because of CS thresholding or 
deletion, some true element gets missed because its value was too 
small, it will, w.h.p., get detected by CS at a future time when I EN 
increases. Also, as long as rank{AT) > size{T) for the currently 
estimated T (which may contain some extra coordinates), the esti- 
mation error will increase beyond MMSE, but will not blow up, i.e. 
it will still converge, but to a higher constant value. 

4. SIMULATION RESULTS 

We simulated a time sequence of sparse ?n=256 length signals, xt, 
with maximum sparsity Smax- Three sets of simulations were run 
with Smax= 8, 16 and 25. The A matrix was simulated as in Q by 
generating n x m i.i.d. Gaussian entries (with n = 72) and normal- 
izing each column of the resulting matrix. Such a matrix has been 
shown to satisfy the UUP at a level Clogm Q. The observation 
noise variance, a^i,^ = {{l/3)^Smax/n)^ (this is taken from Q). 
The prior model on Xt was ^ with — 9 and a^y^ — 1. Ti 

(support set of xi) was obtained by generating Smax — 2 unique in- 
dices uniformly randomly from [1 : m]. We simulated an increase 



in the support at t = 5, i.e. Tt — Ti, < 5, while alt — 5, we 
added two more elements to the support set. Thus, Tt = Ts, > 5 
had size Smax- Only addition to the support was simulated. 

We used the proposed KF-CS algorithm (Algorithm 1) to com- 
pute the causal estimate Xt of xt at each t. The resulting mean 
squared error (MSB) at each t, Ex,y[\\xt — xt\\2], was computed 
by averaging over 100 Monte Carlo simulations of the above model. 
The same matrix. A, was used in all the simulations, but we averaged 
over the joint pdf of J/, i.e. we generated Ti, Ts, {ut)Tt,, wtjt — 
1, ... 10 randomly in each simulation. Our simulation results are 
shown in Fig. 1 (KF-CS is labeled as CS-KF in plots by mistake). 
Our benchmark was the genie-aided KF, i.e. an Smax-ovA&t KF with 
known Ti and Ts, which generates the MMSB estimate of xt- We 
simulated two types of KF-CS methods, one with known Ti, but un- 
known Ts and the other with unknown Ti and T5. Both performed 
almost equally well for Smax ~ 8, but as Smax was increased much 
beyond the UUP level of A, the performance of the unknown Ti 
case degraded more (the CS assumption did not hold). We also show 
comparison with regular CS at each t, which does not use the fact 
that Tt changes slowly (and does not assume known Ti either). This 
had much higher MSB than KF-CS. The MSB become worse for 
larger Smax- We also implemented the full KF for the 256-dim state 
vector. This used ^ with Qt = cr^j^g J256X256, i-e. it assumed no 
knowledge of the sparsity. Since we had only a 72-length obser- 
vation vector, the full system is not observable. Since all non-zero 
modes are unstable, its error blows up. 



5. CONCLUSIONS AND FUTURE DIRECTIONS 



To the best of our knowledge, this is the first work on extending 
the CS idea to causally estimate a time sequence of spatially sparse 
signals. We do this by using CS to estimate the signal support at the 
initial time instant, followed by running a KF for the reduced order 
model, until the innovation or filtering error increases. When it does, 
we estimate the "change in support" by running CS on the filtering 
error. This has much lower error since the "change" is much sparser 
than the actual signal. Open questions to be addressed in future are 
(a) the analysis of the stability of KF-CS, (b) comparison of KF-CS 
error with that of regular CS, (c) studying how and when to delete 
constant coefficients, (d) KF-CS for compressible signal sequences, 
and (e) its extensions to large-dimensional particle filtering 1141 . 
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